#load an R package
library(rdacca.hp)
library(dplyr)
library(ggplot2)

#===PF==============================
mydata <- read.csv("PF_C_positive.csv")
mydata <- na.omit(mydata)
PSF <- mydata$PSF
env <- mydata[,4:11]

#RDA analysis
hp_PSF <- rdacca.hp(mydata['PSF'], env, method = 'RDA', type = 'R2')

# Individual effect
hp_PSF.ie <- data.frame(hp_PSF$Hier.part, check.names = FALSE)

# Total_explained_variation
hp_PSF$Total_explained_variation <- sum(hp_PSF.ie$Individual)
hp_PSF.ie$env <- rownames(hp_PSF.ie)
hp_PSF.ie$variable <- 'PSF'

# Relative effect percentage
hp_PSF.ie$relative_effect <- 100 * hp_PSF.ie$Individual / hp_PSF$Total_explained_variation
combined_ie <- hp_PSF.ie

# Total_explained_variation
combined_ie$Total_explained_variation <- hp_PSF$Total_explained_variation

combined_ie$env <- factor(combined_ie$env, levels = c("Bac_shannon","Bac_niche","Bac_stability","Bac_keystone",
                                                      "Fun_shannon","Fun_niche","Fun_stability","Fun_keystone"))


yanse=c("#FDAE61","#FDD17E","#F3E990","#DDF199","#B3E0A2","#83CDA4","#57B1AB","#3288BD")


p1 <- ggplot(combined_ie, aes(x = relative_effect, y = env, fill = env)) +
  geom_col(width = 0.85) +
  scale_color_manual(values = yanse)+
  scale_fill_manual(values = yanse)+
  theme_bw() +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.ticks.x = element_line(color = 'gray30'),  
        axis.line.y = element_line(color = 'gray30')) +
  labs(y = NULL, x = NULL, fill = '') +
  theme_bw()+
  theme(legend.position = 'none') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  geom_text(data = combined_ie %>% group_by(variable) %>% slice(1),
            aes(y = 1, label = paste("R² = ", round(Total_explained_variation, 2))),
            hjust = 0, vjust = -0.5, fontface = "italic", size = 3, color = "black") +
  xlim(c(0, 65))

p1

dat_select <- mydata[c('Bac_shannon', 'Bac_stability', 'Bac_niche','Bac_keystone',
                       'Fun_shannon', 'Fun_stability', 'Fun_niche','Fun_keystone','PSF')]

dat_select <- na.omit(dat_select)
dat_select <- data.frame(scale(dat_select))

iv <- list(
  Bacteria = dat_select[c('Bac_shannon', 'Bac_stability', 'Bac_niche','Bac_keystone')],
  Fungi = dat_select[c('Fun_shannon', 'Fun_stability', 'Fun_niche','Fun_keystone')] 
)

hp <- rdacca.hp(dat_select['PSF'], iv, method = 'RDA', type = 'R2')
hp

#===AM==============================
mydata <- read.csv("AM_C_positive.csv")
mydata <- na.omit(mydata)

PSF <- mydata$PSF
env <- mydata[,4:7]

hp_PSF <- rdacca.hp(mydata['PSF'], env, method = 'RDA', type = 'R2')

# Individual effect
hp_PSF.ie <- data.frame(hp_PSF$Hier.part, check.names = FALSE)

# Total_explained_variation
hp_PSF$Total_explained_variation <- sum(hp_PSF.ie$Individual)
hp_PSF.ie$env <- rownames(hp_PSF.ie)
hp_PSF.ie$variable <- 'PSF'
hp_PSF.ie$relative_effect <- 100 * hp_PSF.ie$Individual / hp_PSF$Total_explained_variation

combined_ie <- hp_PSF.ie
# Total_explained_variation
combined_ie$Total_explained_variation <- hp_PSF$Total_explained_variation
combined_ie$env <- factor(combined_ie$env, levels = c("AM_shannon","AM_niche",
                                                      "AM_stability","genera"))

yanse=c("#FDAE61","#FDD17E","#F3E990","#DDF199","#B3E0A2","#83CDA4","#57B1AB","#3288BD")


p1 <- ggplot(combined_ie, aes(x = relative_effect, y = env, fill = env)) +
  geom_col(width = 0.85) +
  scale_color_manual(values = yanse)+
  scale_fill_manual(values = yanse)+
  theme_bw() +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.ticks.x = element_line(color = 'gray30'),  
        axis.line.y = element_line(color = 'gray30')) +
  labs(y = NULL, x = NULL, fill = '') +
  theme_bw()+
  theme(legend.position = 'none') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  geom_text(data = combined_ie %>% group_by(variable) %>% slice(1),
            aes(y = 1, label = paste("R² = ", round(Total_explained_variation, 2))),
            hjust = 0, vjust = -0.5, fontface = "italic", size = 3, color = "black") +
  xlim(c(0, 80))


p1




dat_select <- mydata[c('AM_shannon','AM_stability','AM_niche','PSF')]
dat_select <- na.omit(dat_select)  
dat_select <- data.frame(scale(dat_select))


iv <- list(
  shannon = dat_select[c('AM_shannon')],
  stability = dat_select[c('AM_stability')], 
  niche = dat_select[c('AM_niche')]
)

hp <- rdacca.hp(dat_select['PSF'], iv, method = 'RDA', type = 'R2')
hp
